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TITLE OF THE INVENTION 

Automatic Multi-Dimensional Intravascular Ultrasound Image 
Segmentation Metliod 

FIELD OF THE INVENTION 

[0001] The present invention generally relates to image segmentation. 

More specifically but not exclusively, the present invention is concerned with an 
intravascular ultrasound image segmentation technique for characterizing blood 
vessel vascular layers from Intravascular ultrasound image sequences. 

BACKGROUND OF THE INVENTION 

[0002] Over the past few years, intravascular ultrasound (IVUS) 

technology has become very useful for studying atherosclerotic disease. IVUS Is a 
medical imaging technique that produces cross-sectional Images as a catheter is 
pulled-back inside a blood vessel. These images show the lumen but also the 
layered structure of the vascular wall. IVUS provides quantitative assessment of 
the vascular wall, Information about the nature of atherosclerotic lesions as well as 
the plaque shape and size such that In clinic, IVUS was rapidly recognized as a 
valuable tool in diagnosis and in pfe-intervention analysis of atherosclerosis. 

[0003] The ability to characterize the vascular wall was initially 

demonstrated in 1989 by Gussenhoven et al., in "Arterial wall characteristics 
determined by intravascular ultrasound imaging: An in vitro study"' (J. Am. ColL 
Cardiol., vol. 14, no. 4, pp. 947-952, 1989). Also, studies of the mid-90s by IVIintz 
et aL, in "Atherosclerosis in angiographically 'normal' coronary artery reference 
segments: An intravascular ultrasound study with clinical correlations" (J. Am. ColL 
Cardiol., vol. 25, no. 7, pp. 1479-1485, 1995), showed, based on IVUS, that 40% 
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of angiographically normal vessel were in fact atherosclerotic. 

[0004] By using IVUS, it was also demonstrated by Colombo et aL, in 

"Intracoronary stenting without anticoagulation accomplished with intravascular 
ultrasound guidance" {Circulation, vol. 91, pp. 1676-1688, 1995) that conventional 
stent implantation resulted in incomplete apposition and expansion causing 
thrombosis, which had the result of changing the clinical practice. 

[0005] IVUS is also expected to play an . important role in 

atherosclerosis research. For example, as demonstrated by Nissen et al., in 
"Application of intravascular ultrasound to characterize coronary artery disease 
and assess the progreission or regression of atherosclerosis" {Am. J. CardioL, vol. 
89, pp. 24B-31 B, 2002), IVUS helps to achieve precise evaluation of the disease 
in new progression7regression therapies. Experts also generally agree that IVUS 
imaging adds precious complementary information to angiography which only 
shows a projection of the lumen, as taught by Nissen et al., in "Intravascular 
ultrasound : Novel pathophysiological insights and current clinical applications" 
{Circulation, vol. 1 03, pp. 604-61 6, 2001 ). 

[0006] Over the last few years, new signal processing strategies were 

applied to IVUS signals for extracting information on the elastic properties of the 
vascular wall. For example, a new imaging technique named "intravascular or 
endovascular ultrasound elastography" was proposed by de Korte et al., in 
"Intravascular elasticity imaging using ultrasound - Feasibility studies in phantoms" 
{Ultrasound Med, Biol., vol. 23, pp. 735-7^, 1997). Recently, Brusseau et al. in 
"Fully automatic luminal contour segmentation in intracoronary ultrasound imaging 
- A statistical approach" {IEEE Trans. Med. Imag., vol. 23, pp. 554-566, 2004) 
suggested to use a pre-segmentation of the structures of the vascular wall 
identified from IVUS images to help implementing elastography algorithms. This 
constitutes another Important domain of application of IVUS multi-dimensional 
image segmentation. 
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[0007] The tomographic nature of IVUS makes 3D reconstruction of the 

vessel wall possible. Furthemipre, 2D and 3D quantitative measurements of 
atherosclerotic disease such as plaque volume, intima-media thickness, vascular 
remodeling, and lumen area stenosis can be retrieved from IVUS data as 
disclosed by Mintz et al., in "American college of cardiology clinical expert 
consensus document on standards for acquisition, measurement and reporting of 
intravascular ultrasound studies (IVUS)" (J^ Am. poll. CardioL, vol. 37, no. 5, pp. 
1478-1492,2001). 

[0008] However, a typical IVUS acquisition generally contains several 

hundreds of images, which has the effect of making analysis of the data a long and 
fastidious task that is further subject to an important variability between intra- 
observers and inter-observers when non-automatic methods are used. These 
aspects generate important constraints against the clinical use of IVUS. Other 
constraints related to the use of IVUS Include poor quality image due to speckle 
noise, imaging artifacts, and shadowing of parts of the vessel wall by calcifications. 

[0009] So far, a number of segmentation techniques have been 

developed for IVUS data analysis and were introduced to overcome the 
hereinabove discussed drawbacks. Generally speaking, a portion of these 
techniques are based on local properties of image pixels, namely with the gradient- 
based active surfaces as introduced by Klingensmith et al., in "Evaluation of three- 
dimensional segmentation algorithms for the identification of luminal and medial- 
adventitial borders in intravascular ultrasound images" (IEEE Trans. Med. Imag., 
vol. 19, no- 10, pp. 996-1011, 2000) and the pixel intensity combined to gradient 
active contours as introduced by Kovalski et a!., in. "Three-dimensional automatic 
quantitative analysis of intravascular ultrasound Images" (Ultrasound Med. Biol., 
vol. 26, no. 4, pp. 527-537, 2000). 

[0010] Graph search was also investigated using local pixel features, 

for Instance, with Sobel-like edge operator as disclosed by Zhang et al., In "Tissue 
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characterization in intravascular ultrasound images" {IEEE Trans. Med. Imag., vol. 
17, no. 6, pp. 889-899, 1998) and with gradient associated to line patterns 
con-elation as demonstrated by Von Birgelen et al., in "Morphometric analysis in 
three-dimensional intracoronary ultrasound : An in vitro and in vivo study using a 
novel system for the contour detection of lumen and plaque" (Am. Heart J., vol. 
132, no. 2, pp. 516-527, 1996). 

[0011 1 The other portion of the IVUS segmentation work was based on 

more global or region infomfiation. For instance, texture-based morphological 
processing was considered as disclosed by Mojsilovic et al., in "Automatic 
segmentation of Intravascular ultrasound images: A texture-based approach" (Ann. 
Biomed. Eng., vol. 25, no. 6, pp. 1059-1071, 1997). Gray level variances were 
then used for the optimization of a maximum a posteriori (MAP) estimator 
modeling ultrasound speckle and contour geometry as demonstrated by Haas et 
al., in "Segmentation of 3D intravascular ultrasonic Images b^sed on a random 
field model" {Ultrasound Med. Biol., vol. 26, no. 2, pp. 297-306, 2P00). 

[0012] In addition, some studies defining only the lumen boundary and 

not using the full IVUS potential can be found In the literature. Still, in 2001, the 
clinical expert consensus from the American College of Cardiology in the 
hereinabove cited document by Mihtz et al. reported that no IVUS edge detection 
method had found widespread accisptance by clinicians. 

[0013] Recently, graph search was revisited using other edge filters, as 

disclosed by Koning et al., in "Advanced contour detection for three-dimensional 
Intracoronary ultrasound: A validation - in vitro and In vivo" {Int. J. Cardiac Imag., 
vol. 18, pp. 235-248, 2002). 

[0014] Other recent models and methods were proposed, such as 

elliptical template fitting as demonstrated by Welchert et al., in "Virtual 3D IVUS 
model for Intravascular brachytherapy planning: 3D segmentation, reconstruction, 
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and visualization of coronary artery architecture and orientation" {Med. Phys., vol. 
30. no. 9, pp. 2530-2536, 2003) and multiagent segmentation by Bovenkamp et 
al., in "Multiagent IVUS image interpretation" {SPIE Proceedings : Medical Imaging 
2003 : Image Processing, vol. 5032, San-Diego, USA, 2003, pp. 619-630). 
However, tliese new models were again using local pixel or edge information and 
they were not taking advantage of the statistical information of IVUS data (speckle 
texture). 

[0015] Since image pixels in IVUS have pixel gray values generally 

distributed according to a Rayleigh probability density function (PDF) in B-mode 
(brightness modulation) imaging of uniform scattering tissues, as demonstrated by 
Wagner et al., in ''Statistics of speckle in ultrasound B-scans" {IEEE Trans. Son. 
Ultrason., vol. 30, no. 3, pp. 156-163, 1983), it is believed that PDF features can 
be of value for IVUS segmentation. This infonmation is hypothetically more suitable 
for IVUS image analysis, especially when contrast is low between layers of the 
vascular wall. In addition, because the IVUS radio-frequency (RF) mode generally 
provides a better spatial resolution than B-mode imaging, it is also expected that 
the Gaussian PDF of RF images can be exploited for image segmentation. 

[0016] Since the atherosclerotic plaque structure on the vascular wall 

can have an irregular and complex shape that is rarely elliptical, a fast marching 
method as disclosed by Sethian in "Level Set Methods and Fast Marching 
Methods: Evolving Interfaces in Computational Geometry, Fluids Mechanics, 
Computer Vision and Materials Science" (2nd ed. Cambridge, UK: Cambridge 
University press, 1999) and by Osher et al., in "Fronts propagating with curvature- 
dependent speed: Algorithms based on hamilton-jacobi formulations" {J. Comput. 
Phys., vol. 79, pp. 12-49, 1988), can be used to handle topological changes and 
contour irregularities generated by IVUS images. Further, the fact that a fast 
marching method propagates interfaces in the direction of the boundaries through 
an exhaustive analysis of the propagation region has the effect, of decreasing the 
variability of segmentation results. 
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SUMMARY OF THE INVENTION 

[0017] More specifically, in. accordance with the present Invention, 

there is provided an image segmentation method for estimating boundaries of 
layers in a multi-layer body, the method including providing Image data of the 
multi-layer body, the Image data representing a plurality of ^image elements. The 
method further includes detenmlning a plurality of initial interfaces corresponding to 
regions of the image datd to segment, and concurrently propagating the initial 
interfaces corresponding to the regions to segment and thereby estimating the 
boundaries of the layers of the multi-layer body. Propagating the initial interfaces 
Including using a fast marching model based on a probability function describing at 
least one characteristic of the image elements. 

[0018] There is furthermore provided an image segmentation method 

for estimating boundaries of layers in a multi-layer fcfody, the method including 
providing image data of the multi-layer body, the image data representing a 
plurality of image elements. The method further includes determining a plurality of 
initial interfaces corresponding to regions of the image data to segment, and 
concurrently propagating the initial interfaces corresponding to the regions to 
segment the regions and estimate the boundaries of the layers of the multi-layer 
body. Propagating the initial interfaces includes using a fast marching model based 
on a gradient function describing at least one characteristic of the image elements. 

[0019] There Is furthermore provided an image segmentation method 

for estimating boundaries or layers in a pulsating multi-layer blood vessel, the 
method Including: providing IVUS image data of the pulsating multi-layer blood 
vessel, determining initial Interfaces corresponding to the regions of the IVUS 
image data to segment, dividing wall pulsations of the IVUS image data, into a 
discrete number of phases with adjustable pulsation phase labels, assigning the 
pulsation phase labels to 2D IVUS frames of the IVUS Image data, dividing the 
IVUS image data according to the phases and propagating the Initial interfaces 
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according to a fast marching model by simultaneously estimating a mixture of 
probability density functions in the IVUS image data for each of the regions to 
segment and according to each of the phasies. 

[0020] There is furthemiore provided an Image segmentation method 

for estimating boundaries of layers in a multi-layer body, the method including: 
providing Image data of the multi-layer body, the image data representing a 
plurality of image elements. The method further includes determining initial 
Interfaces corresponding to the regions of the image data to segment and 
propagating the initial interfaces according to a fast marching model. Propagating 
the initial interfaces Includes, for each region to segment, simultaneously 
estimating a speed function for the propagation of the initial interfaces based on a 
probability function describing at least one characteristic of the image elements, 
and mapping a time function of the propagating initial Interfaces. 

[0021] There is furthermore provided a data acquisition system for 

segmenting images by estimating boundaries of layers in a multi-layer body, 
including: a catheter including a transducer for providing image data representing a 
plurality or image elements and a data acquisition tool including: a digitizer in 
communication with the transducer for digitizing the image data, a memory for 
receiving and storing the digitized image data, a calculator for estimating, for each 
of the layers, probability functions describing at least one characteristic of the 
image elements, a processor in communication with the memory and the calculator 
for simultaneously estimating the boundaries of the layers of the digitized image 
data by using a fast marching model based on the estimated probability functions. 

[0022] The foregoing and other objects, advantages and features of the 

present invention will become more apparent upon reading of the following non- 
restrictive description of illustrative embodiments thereof, given by way of example 
only with reference to the accompanying drawings. 
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BRIEF DESCRIPTION OF THE DRAWINGS 
[0023] In the appended drawings: 

[0024] Figure 1 is a 2D IVUS frame view representing the various 

layers of a blood vessel in image data which is used in a segmentation method for 
detecting layer boundaries according to an illustrative embodiment of the present 
invention; 

[0025] Figure 2 is another 2D IVUS frame view representing one 

irregulariy shaped layer of the blood vessel shown In Figure 1 ; 

[0026] Figure 3 is a chart view showing the gray levels present in IVUS 

image data and the mixture of the probability density functions detected per layer; 

[0027] Figure 4 is a longitudinal view generated from IVUS image data 

showing an operation of a method for detecting the layer boundaries; 

[0028] Figure 5 is a 2D IVUS frame view intersecting the longitudinal 

view of Figure 4 at point A; 

[0029] Figure 6 Is a flowchart schematically representing simulated 

IVUS image data generated from a plurality of 2D IVUS frames as the one shown 
in Figure 1 ; 

[0030] Figure 7a is a cross-sectional view of a simulated blood vessel 

reconstructed according to a segmented 2D IVUS frame; 

[0031] Figure 7b is a 2D IVUS frame view of the simulated blood vessel 

of Figure 7a and generated by the method shown in Figure 6; 
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[0032] Figure 7c is a 2D IVUS frame view sliowing the segmentation 

results on the 2D IVUS frame view of Figure 7b, with the boundaries detected by 
the 3D fast marching method based on probability density functions; 

[0033] Figure 7d Is a 2D IVUS frame view showing the segmentation 

results on the 2D IVUS frame view of Figure 7b, with the boundaries detected by 
the 3D gradient fast marching method; 

[0034] Figure 8a is another cross-sectional view of a simulated blood 

vessel reconstructed according to a segmented 2D IVUS frame; 

[0035] Figure 8b is a 2D IVUS frame view of the simulated blood vessel 

of Figure 8a and generated by the method shown in Figure 6; 

[0036] Figure 8c is a 2D IVUS frame view showing the segmentation 

results on the 2D IVUS frame view of Figure 8b, with the boundaries detected by. 
the 3D fast marching method based on probability density functions; 

[0037] Figure 8d is a 2D IVUS frame view showing the segmentation 

results on the 2D IVUS frame view of Figure 8b, with the boundaries detected by 
the 3D gradient fast marching method; 

[0038] Figure 9a is a 2D IVUS frame view similar to the one shown in 

Figure 1 ; 

[0039] . Figure 9b is a 2D IVUS frame view showing the segmentation 

results on the 2D IVUS frame view of Figure 9a, with the boundaries detected by 
the 3D fast marching method based on probability density functions; 



[0040] 



Figure 9c is a 2D IVUS frame view showing the segmentation 
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results on the 2D IVUS frame view of Figure 9a, with the boundaries detected by 
the 3D gradient fast marching method; 

[0041] Figure 10a is a 2D IVUS frame view similar to the one shown in 

Figure 1 ; 

[0042] Figure 10b is a 2D IVUS frame view showing the segmentation 

results on the 2D IVUS frame view of Figure 10a, with the boundaries detected by 
the 3D fast marching method based on probability density functions; 

[0043] Figure 10c is a 2D IVUS frame view showing the segmentation 

results on the 2D IVUS frame view of Figure 10a, with the boundaries detected by 
the 3D gradient fast marching method; 

[0044] Figure 11 is a longitudinal view showing a volumic 

reconstruction of the vessel layers detected according to the fast-rnarching method 
based on probability density functions; 

[0045] Figure 12a is a 2D IVUS frame view representing the various 

layers of a blood vessel on simulated RF image data; 

[0046] Figure 12b is a 2D IVUS frame view showing the segmentation 

results on the RF image data of Figure 12a, with the boundaries detected by the 
3D fast-marching method based on probability density functions; 

[0047] Figure 13a is a detailed schematic view showing a first example 

of propagation area for detecting a layer boundary using the fast marching method; 

[0048] Figure 13b is a detailed schematic view showing another 

example of propagation area for detecting a layer boundary using the fast 
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marching method; 

[0049] Figure 14a is a 2D IVUS frame view which is undersampled with 

respect to a typical 2D IVUS frame and which may be used according to a second 
illustrative embodiment of the present invention; 

[0050] Figure 14b is another 2D IVUS frame view which is 

undersampled with respect to a typical 2D IVUS frame but with a higher resolution 
than the 2D IVUS view shown in Figure 14a; 

[0051] Figure 14c is another 2D IVUS frame view which is 

undersampled with respect to a typical 2D IVUS frame but with a higher resolution 
than the 2D IVUS view shown in Figure 14b; 

[0052] Figure 14d is a typical 2D IVUS frame view with a higher 

resolution than the 2D IVUS views shown in Figure 14a to 14c; 

[0053] Figure 15 is a detailed schematic view showing template regions 

of a vessel for detecting layer boundaries in a method according to a third 
illustrative embodiment of the invention; 

[0054] Figure 16 is a cross-sectional longitudinal view of an IVUS 

image data showing a sawtooth artifact typically caused by pulsations of blood 
vessels; 

[0055] Figure 17 Is a chart view showing layer areas resulting from 

segmentation detected by a 3D fast marching method based on probability density 
functions according to a fifth illustrative embodiment of the present invention; 



[0056] 



Figure 18 is a flowchart representing a method for determining 
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the initialization contours according to tlie first illustrative embodiment of the 
present invention; 

[0057] Figure 19 is a block diagram representing an in-vivo IVUS data 

acquisition and processing device according to the first illustrative embodiment of 
the present invention; 

[0058] Figure 20 is a flowchart representing an IVUS method for 

initializing propagating interfaces generated from lower resolution segmentation 
results according to the second illustrative embodiment of the present invention; 

[0059] Figure 21 is a flowchart representing an exploration method 

from a wide propagating area at low resolution to a reduced propagating area at 
high resolution according to the second illustrative embodiment of the present 
invention; 

[0060] Figure 22 is a flowchart representing a template region 

searching method for automatically finding initial interfaces of the layers according 
to the third illustrative embodiment of the present invention; 

[0061] Figure 23 is a flowchart representing an automatic estimation 

method of the probability density function mixture parameters based on the 
iterative conditional estimation according to the fourth illustrative embodiment of , 
the present invention; and 

[0062] Figure 24 is a flowchart representing the method for using layer 

pulsation assessment in the boundary detection segmentation process of the fast- 
marching method according to the fifth illustrative embodiment of the present 
invention. 
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DETAILED DESCRIPTION OF THE ILLUSTRATIVE EMBODIMENTS 

[0063] The non-restrictive illustrative embodiments of the present 

invention relate to a method and device for concurrently estimating boundaries 
between the plurality of layers of a blood vessel from IVUS image data. The 
method and device involve a segmentation of the IVUS image data by propagating 
interfaces in each layer to be estimated from initial interfaces that are generated 
from the IVUS image data. The technique for estimating the boundaries of the 
various layers uses a fast marching method based on probability functions, such 
as for example a probability density function (PDF) or gradient function to estimate 
the distribution color map of images, such as for example to estimate the gray 
levels or the multi-colored levels In images. 

[0064] The following description is organized as follows. First of all, a 

PDF estimation technique for the different vessel layers will be presented. Then, 
an IVUS 3D fast marching method based on the estimated PDFs and based on the 
gray level gradient will be explained and followed by an initializing technique. 
Finally, segmentation results on experimental B-mode data, simulated B-mode and 
simulated RF data and will be reported and discussed. 

[0065] IVUS images are generally provided from an ultrasound 

transducer at the tip of a catheter that is pulled back inside a blood vessel and 
produces a plurality of IVUS 2D frames. A typical IVUS 2D frame is illustrated in 
Figure 1. As illustrated, the 2D frame of Figure 1 shows the catheter and some 
layers of the blood vessel such as, for example, the lumen, the intima and plaque, 
the media and the surrounding tissues. Figure 2 illustrates how the boundary of the 
lumen may be irregularly shaped. 

[0066] The IVUS 2D frames are ultrasonic images made from a 

plurality of pixels generally colored with various shades of gray. In B-mode 
(brightness modulation) or RF-mode (radio-frequency) imaging, such as for 
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example in IVUS data, a Rayleigh or a Gaussian probability density function (PDF) 
can be used, respectively, to model the color map distribution of tlie ultrasonic 
speckle pattern in a uniform scattering tissue. Wlien more than one layer of tissue 
is present, the color map distribution of a whole IVUS Image data can then be 
modeled by a mixture of Rayleigh or Gaussian PDFs, depending on the mode 
selected on tlie instrument. 

[0067] The illustrative embodiment that follows generally considers 

IVUS B-mode imaging, but one ordinary skilled in the art will easily understand that 
similar equations can be provided for Gaussian PDFs if the RF-mode is 
considered. For more details, see IHastie et al., in "The elements of statistical 
learning. Data mining, inference and prediction" (New York, USA: Splnger, pp. 
236-242,2001). 

[0068] In this first illustrative embodiment, a Rayleigh probability 

density function (PDF) Pxi^) models the gray level color map distribution by using 

a parameter , where x is the gray level taking values situated, for example, in 
the range [1, 256]. In this particular example, the Rayleigh probability density 
function (PDF) Is given by equation 1 : 



[00691 a y 2a J (1) 

[0070] with ^»«* >0, and the variance = a^{4-fr)/2 , 

[0071] IVUS data are modeled by a mixture of IVI Rayleigh PDFs 

(corresponding to iVI different layers of the blood vessel) with parameters 

® = {(^7.07)}^=, where is the proportion of the component of the mixture of 
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the M Rayleigh PDFs, so that . The global data PDF mixture then 

becomes : 



[0072] p^3(; 



®)^f,o>jp(x\a)) (2) 



[0073] To describe the PDF mixture for the global IVUS data, the 

parameters (^y^y) of each PDF composing the mixture need to be estimated. In 
IVUS data, the occurring probability of the gray level values x, or observed data, 
can generally be measured by computing the image histogram, as shown in Figure 
3, but the blood vessel layer to which each pixel of an IVUS image belongs is 
generally unknown or hidden for images that are not already segmented. 

[0074] The Expectation-Maximization algorithm (EM) is an iterative 

computation technique of maximum likelihood estimates for incomplete data, as 
presented by Dempster et aL, in "Maximum likelihood from Incomplete data via the 
EM algorithm" (J. Roy. Stat Soc. S, vol. 39, no, 1, pp. 1-38, 1977), which can be 
used to provide the unknown parameters or hidden information of the probability 
density functions (PDFs). Because the IVUS data are incomplete in terms of 
maximum likelihood estimation, the EM algorithm can be applied to evaluate the 
missing or hidden mixture parameters of the Rayleigh or Gaussian PDFs. 



[0075] The EM algorithm therefore helps to describe the global data 

PDF mixture because ©, a mixture parameter maximizing the likelihood of 

p(x\&) , cannot be solved analytically. A hidden variable Y representative of the 

tissue class (vascular layer to which the pixel belongs) and taking values situated 
in the range [1, ...,M], must be introduced at this point. The log-likelihood of the 

joint distribution of (x,Y)={{x^,yg)}f^^, where N represents the size of the IVUS 

data, is: 
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[0076] log(px.Yie(^.:K I ©))= i;iogi>(y,)p(x, I y,,&) (3) 

[0077] The first step of the EM algorithm is called the Expectation Step 

which calculates the cost function .Q{®,&) = Ey\log{p{X,Y\&))\X,&'], the 
expected value of the log-likelihood of (X,Y), the joint distribution, given the 

observed data X and ©' = )}^,. a previous estimate of the PDF mixture 

parameters. 



[0078] The next step is to determine a new estimate © of the PDF 

mixture parameters by maximizing Q{&,@') with respect to parameters © ; this 
operation can now be performed analytically. 

[0079] The detailed PDF parameter estimation procedure via the EM 

algorithm is therefore: 

[0080] • Initialize ©' , the previous estimate of PDF mixture parameters, 

[0081] • Expectation Step: 

[0082] Evaluate the cost function : 

[0083] e(©,eO = iEY[log(p(X, Y 1 ©))| X,©'] (4) 



[0084] 



M N > . . . 

/=1 M 



(6) 
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[0085] 



Calculate pG, = y | x , , G ') = . 

E ^kPk(x, la^ 



*=i 



according to Bayes aile, and using the previous parameter estimate 0' and 
Equation 1 . 



[0086] 



Maximization Step: 



[0087] 

parameters : 



Calculate ©, the new estimate of the PDF mixture 



[0088] 



Q>y =argmax^ I q(@,&)+Z l-^^i 



M 



[0089] 



(6) 



[0090] 

to1. 



where A = N is a Lagrangian making the sum of the(27y iequal 



[0091] dj = argmax^j e(©,©') 



[0092] 



a, = 



(7) 



[0093] • If @^&\ update previous estimate ©' = ©, and repeat the 

Expectation and IVIaximization steps. 
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[0094] In summary, the EM. algorithm maximizes the likelihood of the 

joint distribution of the observed and hidden data by estimating the posterior 
distribution with pYfx..&iy\X'®')- An interesting property of the EM algorithm is that 
it is guaranteed that the likelihood of the observed data X increases at each 
iteration. 

[0095] For computation efficiency, the EM algorithm is generally 

applied to a randomly drawn subset of the observed data X, which are, In this 
case, a portion of the pixels from the whole IVUS data. For instance, the subset 
size may be about 400 000 pixels when a complete IVUS pullback generally 
contains over 80 000 000 pixels. 

[0096] EM algorithms are otherwise well known to those of ordinary 

skill in the art and, accordingly, will not be further described in the present 
specification. 

[0097], The estimated gray level PDFs of the blood vessel layers can 

then be used to establish a segmentation model in the fast marching framework. 
The fast marching method Is derived from the level-set niodel disclosed by Sethian 
in 'Level Set Methods and Fast Marching Methods: Evolving Interfaces in 
Computational Geometry. Fluids Mechanics, Computer Vision and Materials 
Science" (2nd ed. Cambridge, UK: Cambridge University press, 1999) and by 
Osher et al., in "Fronts propagating with curvature-dependent speed: Algorithms 
based on hamllton-jacobi fomnulations" (J. Comput. Phys., vol. 79, pp. 12-49, 
1988). The fast marching method helps to follow interface propagation. 

[0098] In the level-set model approach, an initial interi^ace is defined as 

the zero level of a function ^ of a higher dimension than the interface. The value 

^(x) of a point x = (x1, x2, .... Xn) e9l"is the distance between that point and the 
initiial interface. The function ^ moves in its normal direction according to a spped 
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function F. The evolution of function ^ Is given by the following Equation 8 with 
initial interface ij>{^^t = o) , 

[0099] ' ^ + ir|v^| = 0. (8) 

[0100] The level-set model is applicable to image segmentation by 

interpreting image boundaries as the final position of the propagating interface, as 
disclosed by Sethian In ''Level Set Methods and Fast Marching Methods: Evolving 
Interfaces in Computational Geometry, Fluids Mechanics, Computer Vision and 
Materials Science" (2nd ed. Cambridge, UK: Cambridge University press, 1999) 
and by R. Malladi et ah, in "Shape modeling with front propagation: A level set 
approach" {IEEE Trans. Pattern Anal. Machine Intell., vol. 17, no. 2, pp. 158-175, 
1995). 

[01011 To achieve this, the speed function F is defined in terms of 

image or shape features and should become close to zero when the propagating 
interface meets with Image boundaries. Since the speed value Is near zero, the 
propagating interface stops on the image boundary, which generally ends the 
segmentation process. 

[0102] Fast marching is a particular case of the level-set model. It, 

consists of an Interface evolving under a unidirectional speed function. In this case; 
for a positive speed function, the propagating interface must be inside the region to 
segment (or outside for a negative speed function), because the propagating 
interface does not explore the region located inside the initial interface. 

[0103] In the fast marching formulation, the evolution of the 

propagating interface is expressed in terms of the arrival time T(x) of the interface 
at point (x). The function T satisfies the following Equation 9, stating that the arrival 



wo 2005/048190 



20 



PCT/CA2004/001970 



time difference between two adjacent pixels increases as the velocity of the 
interface decreases. 



the arrival time function T(x). The construction algorithm, as disclosed by Sethian 
in "A fast marching level set method for monotonically advancing fronts." 
{Proceedings of the National Academy of Sciences of the United States of 
America, vol. 93, pp. 1591-1595, 1996), selects the interface point x having the 
smallest arrival time and calculates the arrival times of its neighbors. This is 
repeated until the Interface has propagated across the whole image or until the 
propagating interface Is considered stationary (when the time gradient Is 
sufficiently high). 

. [0106] The level-set and fast marching equations are independent of 

the Interface dimension. On a discrete 3D grid, neighbors' arrival times are 
updated by solving the following approximation of Equation 9: 



[0104] 



|vr|F=r. 



(9) 



[0105] 



The propagation of the interface is done via the construction of 



[0107] 




[0108] 




[0109] 




(10) 



[0110] 



For the x dimension, 
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[0111] 



[0112] 



where A is the grid element size and {i,j\k) is the 3D position 



of the point having its arrival time calculated. Similar definitions apply for 
DfiJj^T and Df;]j.T , in the y and z dimensions. 



and media) must be identified on the IVUS data, image segmentation is 
simultaneously done via a multiple interface extension of the fast marching method 
as disclosed by Sifakis et al., in "Bayesian level sets for image segmentation" (J. 
Vis. Commun. Image R., vol. 13, pp. 44-64, 2002). A speed function is then 
defined for each propagating interface and the T map is built by selecting the point 
with the smallest arrival time value from ail propagating interfaces. 

[0114] Therefore, the fast marching method with multiple propagating 

interfaces enables simultaneous segmentation of different layers of the blood 
vessel. The multiple interfaces directly depict the layered structure of the blood 
vessel and provide that the boundaries do not overlap. 

[011^] In the PDF-based fast marching method, each interface 

associated to a vessel layer evolves at a velocity defined in terms of the PDF P^^^ 
of the corresponding anatomical structure. The propagation speed of the interface 
m^L, where L is the set 1, 2, . . . ,iV^ of the evolving interfaces, is given by 
Equation 1 1 . 



[0113] 



As stated hereinabove, since multiple contours (lumen, Intima 



-1 



[0116] 




(11) 
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[0117] is the gray level value of pixel s at position {fj,k) in image 

^ » Pfni^s) QncI are the measured occurring probabilities of pixel /, in region 

m and /, respectively. Because the occurring probability is more significant for a 
region than for a single pixel, the speed function is calculated over a certain 
number of neighbors, which are pixels located around position {ij,k), such as 
for example, the 26-connected pixels around position {Uj\k). According to 
Equation 11, the m interface velocities will usually be positive and take higher 
values when inside a region having a grayscale distribution close to P„ . 

[0118] As the propagating interfaces approach the boundaries of the 

blood vessel layers, neighbors start to be distributed under other components of 
the PDFs as stated hereinabove, which has the effect of generally increasing 
Pi{l,) and decreasing P„(/J and therefore, decreasing the intert'ace speed. The 

velocity function of Equation 1 1 has a general form that may be used with any 
types of PDF and provides neighborhood averaging. 

[0119] When used for multiple propagating Interfaces, the fast 

marching segmentation method ends when all adjacent propagating interfaces 
meet with their respective boundaries. Propagating interfaces thus evolve until the 
arrival time 7 map is completely built. 

[0120] Since the gray level gradient is a widely accepted image feature, 

comparisons can be made between the hereinabove disclosed PDF 
implementation of the fast marching segmentation and a gray level gradient, 
implementation of the fast marching segmentation. In the latter case, the speed 
function is given by: 



[0121] 



(12) 
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[0122] where is a Gaussian smootliing filter of standard deviation 

a - The speed function of Equation 12 generally propagates interfaces faster on 
low gradient regions. 

[0123] As stated hereinabove, the fast marching segmentation method 

generally requires that the initial Interface is located inside the region to segment. 
This requirement can be achieved with an initialization procedure, in which 
initialization contours are manually traced with respect to data extracted along 
longitudinal planes of the IVUS data, it can also be automatically performed by 
considering a priori information on the whole IVUS data set, as will be further 
described hereinafter. 

[0124] Generally speaking, the step of selecting data along longitudinal 

planes within the IVUS data Is used, instead of using data from a single 2D IVUS 
frame, since longitudinal planes are able to provide information about the whole 
series of data along the length of the blood vessel. Further, the number of 
manually or automatically traced initialization contours on the longitudinal plane is 
Independent of the number of IVUS 2D frames. 

[0125] Initialization contours may be drawn from different numbers of 

longitudinal planes along the blood vessel. As an example, 3 longitudinal planes 
taken at equally spaced angles over 360 degrees may be selected to cut the IVUS 
data volume. The initialization contours provide reference points for generating the 
set of initial interfaces on each IVUS 2D frame, for each of the layers to be 
estimated. This is generally accomplished by attributing respective reference 
points to the IVUS 2D frame corresponding to each initialization contour points. 

[0126] In the illustrative embodiment of Figures 4 and 5, a reference 

point A is taken along to one initialization contour B estimating the lumen boundary 
on the longitudinal plane shown in Figure 4 (Operation 171 of Figure 18). For 
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instance, the reference point A of Figure 4 is transposed on the IVUS 2D frame of 
Figure 5 to help, generate the initial interface of the lunnen on the IVUS 2D frame 
(Operation 172 of Figure 18). IVlore than one reference points generally coming 
from more than one longitudinal plane are then transposed on the IVUS 2D frame 
(Operation 173 of Figure 18). This initialization step is generally done for each 
boundary layer of the blood vessel which needs to be estimated. 

[0127] For each IVUS 2D frame, slightly shrunk splines passing 

through these reference points are computed -and used to generate the initial 
interfaces (Operation 174 of Figure 18). Therefore, using this, procedure, only a 
few longitudinal planes containing IVUS data volume information are required to 
manually generate the Initialization contours in the longitudinal planes and thereby 
to generate the initial interfaces (Operation 175 of Figure 18) required to initialize 
the concurrent segmentation of multiple vessel layers (Operation 176 of Figure 
18), over several hundreds of images which is a typical number for a typical IVUS 
data. 

[0128] The initial longitudinal contours can also alternatively be 

restarted from any previous reference points (Operation 177 of Figure 18). In this 
manner, the user can explore, on-line and easily, sections of an IVUS data that 
were more difficult to interpret on longitudinal planes. 

[0129] Experimental Testings 

[0130] Experimental testings of the hereinabove proposed non- 

restrictive illustrative method of Figure 18 were conducted on a total of 8 in-vivo 
IVUS pullbacks (600 frames/IVUS data) from diseased superficial femoral arteries. 
These experimental testings were performed on 6 different patients before 
undergoing balloon angioplasty. B-mode IVUS data were acquired with a data 
acquisition system, such as for example a Jomed equipment (In-vision gold, 
Helsingborg, Sweden), using a 20 MHz ultrasonic transducer (181 in Figure 19). 
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IVUS 2D. B-mode images of size 10x10 mm were digitized (digitizer 182 in Figure 
19) on 384 x 384 pixel matrices and stored In a memory (183 In Figure 19) using 
the DICOM standard. The acquisition was done at a 10 Images/sec frame rate and 
the catheter pullback velocity was set to 1 mm/sec generating 0,1 mm thick 2D 
slices. Acquisition parameters were set by a clinician to optimize image quality; 
more specifically, the gain varied from 46 to 54 and the grayscale look-up table 
was set to 5. (mage acquisition was not ECG-gated. 

[0131] To evaluate the robustness of the PDF mixture parameter 

estimation of <Dj and aj, the hereinabove described EM algorithm was run 10 

times in a calculator (184 in Figure 19) for one IVUS catheter pullback, with 
different subsets of pixels at each run of the algorithm. Average PDF mixture 
parameters and standard deviations were calculated for the detected Rayleigh 
PDFs. Since pixel subsets were frorh the same IVUS data, PDF mixture 
parameters were expected to generally converge toward the same values. 

[0132] Once this robustness validation was completed, the EM 

algorithm was applied at the beginning of each segmentation, because PDF 
mixture parameters are specific to each IVUS data, as gain and other parameter 
settings are different between each IVUS data, and as echogenicity of each layer 
is variable from one patient to the other. The detected PDF mixtures were 
composed of 4 distributions (lumen, intima, media, and surrounding tissues), but a 
skilled reader will easily understand that the EM algorithm is general and may 
estimate more PDF distributions of heterogeneous vessel layers If required. 

[0133] Testings were conducted on in-vivo blood vessels and 

numerical simulations of blood vessel IVUS data. The in-vivo B-mode IVUS 
images were segmented with 3D multiple interface fast marching using 
automatically detected gray level Rayleigh PDFs and, as a comparison, using the 
gray level gradient. All catheter pullbacks were segmented three times with both 
3D methods using different sets of initial contours. Lumen, Intima (plaque), and 
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media boundaries were obtained. To quantify the variability of boundary detection 
under various initializations, average and Haussdorf point-to-point distances, as 
disclosed by Chalana et a!,, in "A methodology for evaluation of boundary 
detection algorithms on medical images" {IEEE Trans. Med. Imag., vol. 16, no. 5, 
pp. 642-652, 1997), between resulting contours from different initial contour sets 
were calculated. Haussdorf distance represents the worst case since It generates 
the maximum distance between different segmentation results. Average and 
Haussdorf distances directly depict point-to-point contour variations. 

[0134] Detected boundaries from a whole IVUS catheter pullback 

represent a blood vessel In 3D that can be reconstructed. Reconstruction of the 
lumen and media contours was made from a simple, smoothed contour stack (see 
Figure 11). 

[0135] In addition to the above-described in-vivo validation of the 

illustrative embodiment of the segmentation method, numerical simulations of 
IVUS data were conducted to evaluate segmentation accuracy. Since the exact 
geometry of the simulated data is generally entirely known, direct performance 
calculations of the detected boundary with respect to the . exact geometry of the 
simulated data can be obtained. The simulated B-mode IVUS data were first 
segmented using the same procedure as for the in-vivo data, also including 3 
different sets of initial longitudinal view generating initialization contours. Lumen, 
intima (plaque), and media boundaries were obtained. Average and Haussdorf 
point-to-point distances between detected contours and ground truth boundary 
position were calculated for the segraeRtetton results from each set of initial 
contours. 

[0136] Because the simulation method described in Figure 6 allows 

synthesizing both RF and B-mode IVUS data, the hereinabove described fast- 
matching segmentation method was also tested by using automatically detected 
gray level Gaussian PDFs from the 3D simulated RF images. 
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[01371 The image-formation model that was used to simulate IVUS 

data (echograms) is detailed by Maurice et al., in "Adapting the Lagrangian 
speckle model estimator for endovascular elastography: Theory and validation with 
simulated radio-frequency data" (J. Acoust Soc. Am,, vol. 116. pp. 1276-1286. 
2004). Under assumptions such as space-invariance of the imaging system, IVUS 
images were modeled by a convolution operation between the point-spread 
function, which is the equivalent radio-frequency Image of a single ultrasound 
scatterer, and the function describing the acoustic Impedance mismatch of each 
scatterer of the simulated tissue structuriBs composing the IVUS data. In other 
words, the point-spread function expresses the intrinsic characteristics of the 
ultrasound imaging system. > 

[0138] The implementation of the image-formation model was made for 

a 20 MHz transducer with a 60% bandwidth at -3 dB and a beam width of 0.1 mm. 
For the purpose of these simulations, the media was selected 2 times more 
echpgenic than the lumen; the plaque, 1.5 times more echogenic than the media; 
and the surrounding tissues 2 times nriore echogenic than the media. The 
echogenicity can be seen as the image intensity reflecting the acoustic impedance 
mismatch of the scatterers. The signal to noise ratio (SNR) was set to 20 dB. 

[0139] Figure 6 shows the image formation model used to simulate' RF 

and B-mode IVUS data. From real 2D in vivo IVUS images, as the one shown in 
box C, the segmented contours or yessel boundaries (lumen, plaque of the intima, 
media) are created (box D) from manually traced contours on the IVUS 2D image 
in box C. Box E shows the function z(x,y) representing the acoustic impedance 

variations within the 2D frame from box D, and box F shows a function z(r,qp) 
expressing the same acoustic Impedance variations mapped within the 2D frame 
in polar coordinates. Box G shows the polar point spread function h(r, qp) with a 
beam width that increases with depth and element H is a 2D-convolution operator. 
Processing of function z(r,qp) (box F) with the polar point spread function h(r, <p) 
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(box G) through the 2D-convolution operator (element H) produces a simulated 
polar radio-frequency (RF) image l(r,q)) (box J). Box K illustrates the polar B-mode 
image lB(r, cp) that was computed using the Hilbert transform (see element L) of 

l(r,qp), as presented by Kallel et al. in "Speckle motion artifact under tissue 
rotation." {IEEE Trans. Ultrason., Ferroelect, Freq. Contn, vol. 41, pp. 105-122, 
1994). Box M shows the Cartesian B-mode Image or simulated IVUS image Ib(x, y) 

calculated from the polar B-mode image .lB(r, qp). This simulation strategy was 
repeated for the whole Image data of an IVUS catheter pullback within a diseased 
superficial femoral artery. 

[0140] Results and Discussions 

[0141] As stated hereinabove, the EM algorithm was applied 10 times 

on 1 IVUS catheter pullback to evaluate the robustness of the PDF mixture 
parameter estimation. At each run, PDF parameters were estimated on different 
pixel subsets of the same IVUS data (subsets contained approximately 400 000 
pixels). Average mixture parameters for each detected Rayleigh PDF are shown In 
the follo\Anng Table I. An example of automatically detected Rayleigh PDF mixture 
is also shown In Figure 3 in dotted lines for each layer components, with the 
hereinbefore presented IVUS pullback gray level histogram of the whole data set. 

Table I 



Laver 
Comoonent 


0) (%) 




Lumen 


29.40±0.10 


0.6456±0.0021 


Intima and 
Plaque 


20.g6±0.50 


347.70±13.02 



wo 2005/048190 



29 



PCT/CA2004/001970 



Media 


13.55±0.14 


22.68±0.53 


Surrounding 
Tissues 


36.09±0.67 


2294.58±34.01 



[0142] Table I shows that small variations were found between different 

runs of the EM algorithm. It can be stated that mixture detection of the various 
boundary layers Is a robust and stable process with standard deviations of o) and 
af ranging from 0.3% to 3.7% for several runs of the EM algorithm applied on 
different pixel subsets of the IVUS catheter puliback. The EM algorithm was thus 
applied to the 8 available IVUS catheter puilbacl<s to study PDF variability between 
different patients. The results are shown In the following Table II. 

Table II 



Laver 
Comoonent 




a^ 


Lumen 


18.82±10.44 


5.52±12.50 


Intima and 
Plaque 


27.81 ±14.54 


1052.40±1 405.97 


Media 


15.87±3.61 


339.46±817.80 


Surrounding 
Tissues 


37.50±13.82 


2580.49±654.49 



[0143] Because of instrument settings and echogenicities specific to 

the different vascular structures for a given patient, Table II emphasizes the 
generally high variability between mixture parameters of distinct IVUS catheter 
pullbacics. These results suggest that the EM algorithm is capable of fitting various 
Rayleigh PDF mixtures from different patients. 
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[0144] The numerically simulated and in-vivo IVUS images can then be 

segmented with 3D multiple interface fast marching methods using automatically 
detected gray level PDFs and gray level gradient for comparison purposes. For the 
experimental testing, all IVUS catheter pullbacks were segmented three times with 
both 3D methods using different sets of initial interfaces obtained from the 
Initialization contours generated from the longitudinal planes. 

[0145] The results obtained for the simulated segmentation of the IVUS 

images with the detected gray level Rayleigh PDFs method and with the gray level 
gradient method are shown In Figures 7a to 8d. Figures 7a to 7d are concemed 
with a first blood vessel geometry, shown in Figure 7a, and Figures 8a to 8d are 
concerned with a second blood vessel geometry shown In Figure 8a. 

[0146] Figures 7b and 8b respectively show the simulated IVUS cross- 

sectional 2D B-mode images of the first and second examples of a simulated blood 
vessel, which can be obtained from the above-described method schematically 
illustrated in Figure 6. In Figures 7c and 8c, the detected boundaries of the lumen, 
the thickened intima and media are estimated with the detected gray level PDFs 
method and, in Figures 7d and 8d, with the gray level gradient method. 

[0147] The typical simulated IVUS segmentation results shown in 

Figures 7a to 8d illustrate that detected boundaries were very close to the blood 
vessel layer boundaries. They also reveal that the extemal boundary of the media 
Is smoother with the PDF fast marching than the gradient-based method, but that 
the lumen, which can have a rougher surtace, was detected with sufficient details. 
Gradient methods seemed to trace speckle contours on object boundaries, 
because speckle generally has high gray level intensity differences. 

[0148] The following Table III includes the results of the average 

distance (AD) and Haussdorf distance (HD), which is the maximum distance 
between the estimated boundary and the true boundary in mm, between detected 
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boundaries obtained from different initialization steps and ground truth values (true 
boundaries) obtained from the simulated geometry: In this table, FMM refers to the 
fast marching segmentation method. Symbol * indicates a statistically significant 
better performance with p < 0.05 on paired t-test, whereas symbol § refers to a 
statistical significance of p < 0.01 , 



Table III 



Seamentation 


Lumen 


Plaaue 


Media 


Method 


AD (mm) 


HD (mm), 


AD (mm) 


HD (mm) 


AD (mm) 


HD (mm) 


FMM-3D 


0.072 


0.226 


0.061 


0.154 


0.063 


0.164 


PDFs 


±0.062 


±0.074 


±0.038 


±0.046* 


±0.038 


±0.048§ 


FMM-3D 


0.069 


0.197 


0.060 


0.173 


0.063 


0.180 


Gradient 


±0.056 § 


±0.085 § 


±0.044 


±0.050 


±0.044 


±0.052 



[0149] The average and Haussdorf distance were chosen as 

comparison metrics instead of area or perimeter differences because they directly 
depict point-to-point contour variations. As can be seen in Table III, very low 
average and Haussdorf distances values were obtained, for both PDF- and 
gradient-based three-dimensional (3D) fast marching, demonstrating that the 
method is very powerful for simulated B-mode IVUS segmentation. In fact, average 
deviation ranged from 0.060 to 0.072 mm and worst point-to-point distances were 
between 0.154 and 0.226 mm, which is highly satisfactory. Lower Haussdorf 
distances were obtained on lumen boundary with the gradient method (p < 0.01 ) 
because the blood and Intlma interface generally produces bright echoes for which 
the gradient infomnatlon Is significant. However, on less contrasting boundaries 
such as the intlma (plaque) and media interfaces, statistically significant lower 
Haussdorf distances (p < 0.05) were achieved with the PDF-based method. 

[0150] Examples of results obtained with the gray level Rayleigh PDF 

method and with the gray level gradient method for the in-vivo IVUS data are 
displayed in Figures 9a to 1 0c. The lumen, intlma and media detected boundaries 
are presented for a first cross-sectional IVUS Image (Figures 9a) and a second 
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different cross-seetional IVUS image (Figure 10a). 

[0151] In Figures Sb and 10b, the detected boundaries of the lumen, 

the intima and media are estimated with the gray level PDF based fast-marching 
method and, in Figures 9c and . 10c, with the gray level gradient based fast- 
marching method. 

[0152] A qualitative analysis of the gray level PDF and gray level 

gradient fast marching segmentation methods reveals that the detected 
boundaries are very dose to all vessel layers. More specifically, Figures 9a to 1 0c 
show that vessel layer boundaries of In-vivo IVUS images can be identified even if 
the contrast is very low, as illustrated at 4 o'clock for the collateral vessel in Figure 
9a. Also, detected boundaries on Figure 10b and 10c demonstrate that non- 
circular lumen may be detected with fast marching methods. 

[0153] The following Table IV shows the average distance (AD) and the 

Haussdorf distance (HD) between detected boundaries on in-vlyo data for the gray 
level PDF and gray level gradient fast marching methods for different manual 
initializations of the interfaces. 



Table IV 



Seamentation 


Lumen 


Plaaue 


Media 


Method 


AD (mm) 


HD (mm) 


AD (mm) 


HD (mm) 


AD (mm) 


HD (mm) 


FMM-3D 


0.092 


0.270 


0.092 


0.256 


0.092 


0.256 


PDFs 


±0.095 


±0.147* 


±0.078 


±0.102 § 


±0.083 


±0.113* 


FMM-3D 


0.092 


0.317 


0.090 


0.287 


0.085 


0.302 


Gradient 


±0.104 


±0.148 


±0.080 * 


±0.092 


±0.088 * 


±0.107 



[0154] In Table IV, FMM refers to the fast-marching segmentation 

method, symbol * indicates a statistically significant better performance with p < 
0.05 on paired t-test, whereas symbol § means a statistical significance of p < 
0.01. 



wo 2005/048190 



33 



PCT/CA2004/001970 



[0155] Results indicate that gray level PDF fast marching has the 

smallest Haussdorf distances (p < 0.01). which remains under 0.270 mm for all 
boundaries compared to a value of up to 0.317 mm for the gray level gradient 
implementation. PDF fast marching also has relatively small average distances 
between contours, of 0.092 mm and lower, but are significantly higher than intima 
and media average distances obtained with the gray level gradient method (p < 
0.05). However, the differences between these distances are generally small 
(lower than the pixel size). Thus, 3D fast marching detected boundaries have small 
variations when initialized differently and the maximum distance to the closest 
point, representing the worst case, generally stays low. This tends to Indicate that 
the segmentation performance Is good in regions lacking information, for example 
when the boundaries to be detected are covered with catheter ring-down artifacts 
of lost behind calcium deposits. 

[0156] Figure 11 shows a 3D reconstruction of the lumen and media 

contours obtained with gray level PDF 3D fast marching segmentation for which a 
double stenosis In that patient is clearly seen. In the figure, the light gray 
corresponds to the inner portion of the media, whereas the dark gray is the vessel 
lumen. The gray ^ level gradient fast marching method provided similar qualitative 
results (data not shown). 

[0157] As mentioned hereinabove, the PDF-based fast-marching 

segmentation method can further exploit RF data in place of B-mode data. On RF 
data, the EM algorithm generally searches for a mixture of Gaussian PDFs 
describing the different layered structures of the vessel wall on IVUS Images. 

[0158] Figure 12a shows a simulated 2D RF image taken from the 

whole 3D data set, whereas Figure 12b presents an example of segmentation 
results obtained with the PDF method applied on RF IVUS data. Qualitatively, 
similar performance can be appreciated when comparing those results to Figure 
8c. However, quantitatively, better accuracy is generally expected because of the 
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higher resolution of RF images when compared to B-mode data. 

[0159] In a preliminary version of the IVUS fast marching segmentation 

method disclosed by Roy Cardinal et al., in "Intravascular ultrasound image 
segmentation: A fast marching method" {Lecture Notes in Computer Sciences. 
Proceedings of MICCAI 2003: Medical Image Computing & Computer Assisted 
intervention, vol. 2879, 2003, pp, 432-439), a 2D version of the fast marching 
method was implemented. 

[0160] Generally speaking, a 2D IVUS algorithm uses segmentation 

from previous IVUS images of the catheter pullback to correct initial interfaces. The 
2D segmentation model disclosed by Roy Cardinal in the above-mentioned study 
was applied to a small IVUS catheter pullback of 200 images. Depending on the 
IVUS application, any dimensions can be considered for implementing the fast- 
marching PDF or gradient based method. The present multi-dimensional method is 
general and conceptually considers 1D to ND dimensions, where N is the order of 
the method. Note that N = 4 considers time varying 3D IVUS data. 

[0161] Since a bigger IVUS B-mode clinical database was available for 

the present study, the 2D version of the fast marching segmentation was applied to 
all available catheter pullbacks. The 2D implementation of fast marching arrival 
time (from Equation 10) and speed functions (from Equations 11 and 12) is 
generally straightforward. In 2D, 8-connected pixels (26 connected pixels for the 
hereinabove presented method) were used for averaging neighbors in the speed 
function calculation. 

[0162] The following Table V shows the average distance (AD) and the 

Haussdorf distance (HD) between boundaries of the detected vessel layers, from 
different initializations with 2D fast-marching segmentation. As for the 3D fast 
marching method, the results are for automatically detected PDF- and gradient- 
based algorithms. 
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Table V 



Seamentation 
Method 


Lumen 
AD (nam) HD (mm) 


Plaaue 
AD (mm) HD (mm) 


Media 
AD (mm) HD (mm) 


FMM-2D 
PDFs 


0.093 
± 0.096 


0.279 
± 0.149 


0.093 
± 0.078 


0.262 
±0.102 


0.091 
± 0.083 


0.262 
±0.113 


FMM-2D 
Gradient 


0.096 
±0.106 


0.316 
±0.147 


0.095 
± 0.088 


0.299 
±0.100 


0.085 
± 0.090 


0.304 
±0.104 



[0163] A two way analysis of variance with pairwise multiple 

comparisons using the Turkey test was performed on average and Haussdorf 
distances for 2D (Table V) and 3D (Table IV) fast marching. 

[0164] Statistical tests showed that average distances from the 2D fast 

marching detected boundaries in Table V were not different from the 3D fast 
marching results for all blood vessel layers. It should be noted that 2D algorithms 
used segmentation from previous images of the catheter pullback to correct Initial 
contours, which Increased boundary precision. Thus, alternatively, the 
segmentation performance can be increased by combining this type of correction 
strategy In the 3D fast marching method, by using a multi-scale segmentation 
approach to Initialize a higher resolution data set with low resolution segmentation 
results of the same catheter pullback. As for the 3D fast marching method, the 
gray level PDF fast marching In 2D had lower Haussdorf distances than the gray 
level gradient method (p < 0.05). Since good average distance performance was 
achieved with the gray level gradient method in both 2D and 3D fast marching, this 
Infomnation can be added with advantage to the PDF speed function of Equation 
11. 

[0165] A second non-restrictive illustrative embodiment of the method 

and device according to the present invention will now be described. For the sake 
of brevity, only the differences between the method and device according to the 
first non-restrictive illustrative embodiment and the second non-restrictive 
illustrative embodiment are described hereinbelow. 
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[0166] In the second non-restrictive illustrative embodiment, the fast 

marching method has been modified to enhance the efficiency in computing time. 
In the fast marching method in accordance with the first illustrative embodiment, 
several Interfaces simultaneously propagate across the IVUS image data. In the 

propagation process, the propagating Interfaces and their neighboring areas are 

I 

explored in a significantly detailed manner. Generally, in the segmentation method 
as described in the first Illustrative embodiment, all pixels are analyzed and the 
propagation process takes Into account all preceding interface neighbors through 
the arrival time map construction. Computation time Is thus Increased as the Initial 
interface propagates In a larger Initial segmentation area. 

[0167] In the first Illustrative embodiment, the position of the Initial 

interfaces Is calculated from a shrunk version of manually or automatic traced 
initialization contours taken along a longitudinal plane. Examples of the position of 
Initial Interfaces calculated from initialization contours are shown in Figures 13a 
and 13b. 

[0168] The black region represents the unexplored propagating area 

130, the gray pixels on each side of the propagating area 130 correspond to the 
propagating interfaces 132,134 and the arrows 136a, 136b, 136c and 136d 
represent the propagation direction of the propagating interfaces 132,134. The 
dashed line 138 represents the desired boundary to reach and the solid line is an 
example of initial interface 140 from which the initial propagating Interfaces 132, 
134 were calculated. In Figure 13a, the interfaces will not detect the boundary 138 
because the propagation area 130 wa^^ c i&ir^iA wide enough to completely Include 
the boundary 138, as in Figure 13b. However, propagation In Figure 13a will be 
completed faster than the case of Figure 13b. 

[0169] Therefore, to decrease computational load, shrinking can be 

diminished to create a smaller propagating area 130 (Figure 13a). However, 
because the fast marching method propagates an interface 132,134 under a 
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unidirectional speed function (see arrows 136a, 136b or 136c,136d) the boundary 
138 to be detected must be located Inside the propagating area 130 that will be 
explored during the propagation. 

[0170] A compromise between the dimension of the propagating area 

130 and the computation time is sought. With known 2D fast marching 
segmentation method, this problem is generally solved by using segmentation 
results from previous 2D images of the catheter pullback to connect initial 
interfaces. The initial interfaces 140 are then more precise and the propagating 
area 130 can be set with smaller dimensions. 

[0171] In the 3D fast marching segmentation method of the first non- 

restrictive illustrative embodiment, a correction similar to this 2D correction 
principle can be made through a multi-resolution or multi-scale representation and 
segmentation of the IVUS data. An example of such multi-resolution images in 
IVUS data is shown in Figures 14a-14d, where Mower resolution images are 
obtained by undersampling the original IVUS image by 2' (Operation 191 of Figure 
20) where /= 3, 2, 1, 0 is the resolution level corresponding to Fig. 14a, 14b, 14c 
and 14d, respectively. High scale structures, such as for example the lumen, are 
generally emphasized on lower resolution images. 

[0172] The segmentation results of a lower resolution representation of 

the IVUS data are mapped Into the next level of resolution (Operation 192 of 
Figure 20). Boukerroul et al., in "Segmentation of ultrasound Images - 
multiresolution 2D and 3D algorithm based on global and local statistics" (Pattem 
Recognition Letters, 24:779-790, 2003) and Mignotte et al., In "A multlscale 
optimization approach for the dynamic contour-based boundary detection issue" 
(Computerized Medical Imaging and Graphics, 25(3):265-275, 2001) propose 
related concepts. 



[0173] 



These segmentation-mapped results are used to initialize the 
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interface propagation at this higher resolution level (Operation 193 of Figure 20). 
At a low-resolution level, a fast coarse exploration of a wide propagating area is 
performed to bring the propagating interfaces 132,134 closer to the desired 
boundaries 138 (Operation 201 of Figure 21). The propagation area 130 can then 
be reduced at each higher resolution level since the interfaces 132,134 are 
iteratively corrected (Operation 202 of Figure 21). A larger propagation area 130 
can thus be explored in less calculation over all possible resolutions. 

[0174] At a resolution level /, a pixel represents a 2' x 2' block of 

pixels from the driginal resolution image. In Figures 14a to 14d, the undersampling 
that was used resulted in loss of information at low resolution levels. To overcome 
this loss of Information, a multiscale PDF-based velocity function (similar to 
Equation 11) is developed and used (Operation 203 of Figure 21), where P(ls) is 
replaced with the likelihood of the 2' x 2' block of pixels corresponding to U when 
the propagation is done at level /, as given by equation 13 : 

[0175] P(7J=nP(/,,) (13) 

[0176] where bt is the block of 2' x 2' pixels and P(I^^ ) is the occurring 

probability of the gray level value of pixel S/ in the zero resolution image / . 

[0177] The multlresolutlon and multiscale fast marching segmentation 

methods of Figures 20 and 21 generally allows to iteratively improve the accuracy 
of the detected boundaries without increasing the computation time. 

[0178] A third non-restrictive illustrative embodiment of the method and 

device according to the present invention will now be described. For the sake of 
brevity, only the differences between the third non-restrictive illustrative 
embodiment and the first non-restrictive illustrative embodiment will be described 
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hereinbelow. 

[0179] In this third non-restrictive iliustrative embodiment, the fast 

marching method has been modified to automatically find the Initial interfaces for 
the layers of the vessel wall (lumen, Inside and outside contours of the media) by 
using likelihood maps of. each components of the vessel wail (lumen, intima and 
plaque, media, and surrounding tissues), which are calculated according to the 
estimated PDF mixture parameters. This approach can be seen as an alternative 
to the manual initializations of the vessel Interfaces described hereinabove, 

[0180] The initialization procedure generally finds a rough estimate of 

the true boundaries of the layers that will be further refined into accurate wall 
contours with the multiresolution or multiscale gray level PDF fast marching 
method - 

[0181] The initialization procedure generally starts on a small subset of 

contiguous 2D IVUS frames from the whole catheter pullback (Operation 21 1 of 
Figure 22), For instance, the subset may contain Nfnn 2D IVUS frames to get as 
much infomiation as possible (Ni„it = ten (10) for example), while keeping enough 
correlation between the 2D IVUS frames. For better results, the selected 2D IVUS 
frames are generally of good quality, with no calcification shadows and with a 
generally homogeneous plaque, In order to maximize the available Information for 
detemnining the initial Interface, 

[0182] To find these good quality 2D IVUS frames, a degree of fitting Is 

first calculated between each of the individual frame histogram and the pullback 
PDF mixture (see Figure 3). The degree of fitting can be measured by the 
Kolmogorov-Smirnov test that calculates the distance between the PDF mixture 
and the normalized histogram of an IVUS 2D frame. This test is used to determine 
if an empirical distribution (IVUS histogram) differs from a theoretical distribution 
(mixture). The Kolmogorov-Smirnov distance, between the subset init of Ninit 



wo 2005/048190 



40 



PCT/CA2004/001970 



images and the global IVUS data PDF P. that should be minimized is : 

[0183] ZZI^4r^-^(^)l 04) 

[0184] where U is the gray level value of pixel s in image /; kjO^) is the 

number of pixels having the gray level value U in the J image of the subset; and N 
is the number of pixels in image /. 

[0185] The contiguous frames having the smallest Kolmogorov- 

Smirnov distances are generally chosen. Since the mixture is calculated over the 
whole catheter pullback, it represents the average lumen and blood vessel wall. 
This test thus generally selects the frames that are similar to the average catheter 
pullback, and these frames are used to start the calculation of the initial interface. 

[0186] Calculation of the initial interface is initiated with the search of 

an inner generally elliptical ring shaped region corresponding to the media 
structure 150 (Operation 212 of Figure 22), as shown in Figure 15. This media 
region 150 is distributed according to the hereinabove disclosed media PDF and 
enclosed between a surrounding tissue region 152 and a plaque region 154. The 
media region 150, fhe surrounding tissue region 152 and the plaque region 154 
are generally of fixed size to simplify the initialization procedure and because only 
a rough estimate of the media is generally necessary. The surrounding tissue 
region 152 and the plaque region 154 do not have to represent the whole plaque 
and surrounding tissues since they are defined to provide additional information 
about the media region 150. 

[0187] The initialization procedure generally begins with the search of 

the media region 150 of the blood vessel because .it is believed that the elliptical 
constraints .are more easily assumed for this layer. It was indeed reported by 



wo 2005/048190 



41 



PCT/CA2004/001970 



Hagenaars et al., in "Gamma radiation induces positive vascular remodeling after 
balloon angioplasty: A prospective, randomized intravascular ultrasound scan 
study" (Journal of Vascular Surgery, 36(2):31 8-324, 2002) in which 15 patients out 
of 16 had a dissection after angioplasty of the femoropopliteal artery making the 
lumen Irregularly shaped. 

[0188] Since IVUS data acquisition is often conducted in 

atherosclerosis treatment trial in which patients undergo angioplasty, the 
Irregularities of the layers should be taken into account. Also, the search of the 
initial lumen interface Is reduced to the inside region of the media (Operation 213 
of Figure 22), which generally prevent the propagating interface from leaking into 
collateral branches when they are present. Moreover, the elliptical shape of the 
Initial media region 150 generally produces a closed initial Interface, even If the 
media is hidden behind a calcification shadow. 

[0189] In order to find the media region 150 in the subset of initial IVUS 

frames, an energy function must be associated with the template 158 of Figure 15. 
In the deformable model framework disclosed by Jain et al., in "Deformable 
template models: A review" (Signal Processing, 71(2):109-129, 1998) and by 
Zhong et al., in "Object tracking using deformable templates" (In Sixth International 
Conference on Computer Vision; pages 410-445, 1998), the energy function is 
generally defined to be minimal when the template fits the searched region. 
Deformations are applied to the template 158 to achieve a minimum energy 
function on the image. 

[0190] The energy function that should be minimized to find the media 

region 150 is given by the following Equation 15: 

1 Ninit 

. [0191] ^(^.^)=-E(irS Ziog^(^/k)) (15) 
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[0192] Where Pf// \n) is the occurring probability of the gray level value 

of pixel s in they ^ IVUS image of the initial subset according to the PDF of region 
n; R = {rp,r„,r,\ are the plaque region (r^) 154, media region (r„) 150 and tissue 
region (r, ) 152 of the template 158; Ninn is the size of the initial frame subset. 

[0193] Generally only linear transformations such as for example 

translations, stretchings, and rotations are applied to the template 158 because 
only a generally rough estimate of the media region 150 is needed. The template 
158 fitting is performed at a reduced resolution level / = 1, as described in the 
second non-restrictive illustrative embodiment, in order to minimize the 
computation time, while keeping a large enough media to work with. Different 
known minimization algorithm can be used for minimizing the deformation model. 

[0194] The lumen region (not shown in Figure 15) is then searched 

from the defined initial media boundary. A region without geometrical shape 
constraints is computed: the lumen likelihood map is calculated and a lumen 
region is grown or propagated using this map (Operation 214 of Figure 22). The 
lumen region generally starts at the IVUS catheter that Is inside the lumen and 
generally located in the vicinity of the center of all pullback frames. 

[0195] The lumen region grows by adding the pixels that are most likely 

to be inside the lumen according to the occurring probability, for example if the log- 
likelihood - ^logP(// I r/) of pixel s in the Nina image subset according to the PDF 

of the lumen region n is low enough. The region is generally forbidden to grow 
beyond the boundary of the media region 150. 

[0196] The media and lumen regions are then adjusted or fitted to the 

Ninit contiguous 2D IVUS frames that were used in the initial subset (Operation 215 
of Figure 22). Linear transformations may be applied to these media and lumen 
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regions to maximize their likelihood to each of the contiguous 2D IVUS frames. 
This step is usually performed as for the hereinabove described media defonnable 
template, but the initial intlma and plaque region 154 is bounded by the lumen 
region and the media region 150, and the initial surrounding tissues region 152 is 
bounded on only one side by the media region 150. 

[0197] This procedure is generally repeated for the next subset of 

contiguous 2D IVUS frames In the catheter pullback. However, the process for 
each contiguous 2D IVUS frame generally starts vAth the results of the previous 
defined media template 158. The growth of the lumen region generally starts from 
a shrunk version of the previous average lumen region. The whole IVUS Image 
pullback is therefore initialized in that manner. 

[0198] Alternatively, the segmentation fast marching method as 

described in the first illustrative embodiment may use a combination of the gray 
level gradient information and the gray level PDF information in the calculation of 
the initial interfaces if using only the gray level PDF information turns out to be 
insufficient to generate an automatic initialization as described in the third 
illustrative embodiment. The gray level gradient information could also be 
Integrated to the interface velocity function of Equation 1 1 . 

[0199] In the case of very low-quality images or high ultrasound 

attenuation limiting penetration within the vascular wall, the proposed initial 
boundary calculation procedure might fall to find some initial contours or region 
boundaries. For these particular cases, minimal user interaction might be required 
to conrect some regions of the interfaces. If necessary, this interaction may further 
be included in the segmentation process to minimize the occurrences of having to 
re-segment accurately found boundaries. 

[0200] In the case where a single boundary is available, such as for 

example the boundary between the lumen and the intima, the elliptical template 
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may be modified in such a way as to provide a two (2) region template generally 
corresponding to the regions between the single boundary. The energy function of 
equation 15 thus becomes generally restrained to two (2) regions and the 
remaining initialization procedure remains generally similar to the case of the 
multiple boundary Initialization. 

/ 

[0201] A fourth non-restrictive illustrative embodiment of the method 

and device according to the present Invention, will now be described. For the sake 
of brevity, only the differences between this fourth Illustrative embodiment and the 
first illustrative embodiments will be described. 

[0202] In this fourth non-restrictive illustrative embodiment, the fast 

marching method has been modified to replace the EM local algorithm presented 
in the first illustrative embodiment. The EM algorithm is a local algorithm in which 
the. neighbors information is missing. This information is generally required, such 
as for example, in the case of heterogeneous plaque where the PDFs are 
generally more difficult to estimate. In addition, convergence is generally very slow 
with the EM algorithm such that it can take a significant number of iterations in 
order to be able to estimate the mixture parameter © of Equation (2). 

[0203] The automatic Initial contour detection procedure presented In 

the previously presented third Illustrative is based on the IVUS PDF Information. In 
cases where the EM algorithm cannot be used, the iterative conditional estimation 
(ICE) algorithm that was previously proposed for the mixture parameter estimation 
of incomplete data by Pleczynski In "Champs de markov cach§s et estimation 
conditionnelle Iterative" (Traitement du Signal, 11(2):141-153, 1994) generally 
represents a more robust algorithm, which generally converges faster than the EM 
algorithm. 

[0204] In the PDF mixture estimation presented in the first illustrative 

embodiment, the random variables (X, Y) are referred to as the complete data 
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Where Y is the gray level taking values in [1, .... 256] (observed data), and X is the 
tissue label talcing values [1, M] for a PDF mixture modeling M different tissues 
(hidden information). For the set of pixels S, the realization j' = (y,)„s Y the 
IVUS B-mode or RF Image and x = {x,\^ are the unknown pixel labels. The EM 
algorithm is considered local because the labels x, are considered Independent. In 

the ICE algorithm, X is supposed Markovian i.e. Px (x) is defined with respect to 
the following neighborhood energy function: 



[0205] 



P{x)= exp 



\ ('•') J 



(16) 



[0206] where ^ is an energy function and the summation is for all pairs 

of pixel neighbors {s,t). 

[0207] The first operation of the ICE algorithm is to simulate, such as 

for example with the Gibbs sampler, n realizations (jc,,...,x„) of X according to the 

posterior distribution i^|K.a{-«|;'.®'). with ©' the initial or the previous iteration 

estimate of the PDF mixture parameter © (Operation 221 of Figure 23). The 
posterior distribution is computed using Bayes .rule from the known Px(x) and 

PypceCyl^. ®')'> and the neighborhood Vs of pixel s : 



[0208] 



exp -Yi^ix„x,) -log p^Xys 

^ J 

2 exp -^4>i.x, =k,x,)-\ozp^,{y^ \x,=k) 

*-l teV, 



(17) 



[0209] 



With these simulations of the hidden data, n sets 
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((xi,j;),...,(3(:„,>^)) of complete data are available. 

[0210] The next operation (Operation 222 of Figure 23) is to calculate 

the new value of 0'=— [®(a:i,j/)+...+©(jc„,>')], where © is a parameter estimator 

of the complete data (maximum likelihood for example). Operations 221 and 222 of 
Figure 23 are generally repeated until convergence of the mixture parameter 
estimate is achieved (Operation 223 of Figure 23). 

[0211] For the Rayleigh mixture, It is assumed that each layer structure 

of the B-mode images is a generally uniform scattering tissue with a significantly 
large number of diffusers because the Rayleigh PDFs model the gray level 
distribution of the ultrasound signal under that condition, as disclosed by Wagner 
et al., in "Statistics of specl^le in ultrasound B-scans" (IEEE Transactions on Sonics 
and Ultrasonics, 30(3): 156-1 63. 1983). The same reasoning applies to Gaussian 
PDFs describing RF IVUS images. 

[0212] In the case of highly heterogeneous plaque layer of a diseased 

patient, the Rayleigh or Gaussian PDF might not be sufficient to model the pixel 
gray level distribution. Distributions other than Rayleigh or Gaussian have been 
investigated in modeling of the ultrasound B-mode envelopes or RF signals, 
respectively: Rician distribution as disclosed by Wear et al., in "Statistical 
properties of estimates of signal-to-noise ratio and number of scatterers per 
resolution cell" (Joumal of the Acoustical Society of America, 102(1):635641, 
1997), K distribution as disclosed by Dutt et al., in "Statistics of the log-compressed 
echo envelope" (Journal of the Acoustical Society of America, g9(6):381 7-3825; 
1996) and Nakagami distribution as disclosed by Shankar In "A general statistical 
model for ultrasonic backscattering from tissues" (IEEE Transactions on 
Ultrasonics, Ferroelectrics, and Frequency Control, 47(3):727- 736, 2000). 



[0213] 



The ICE algorithm generally has no limitation for the type of 
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statistical distribution to be modeled, as long as a parameter estimate of the 
complete data can be calculated. Moreover, if a model of mixed distribution types 
is necessary, the generalized ICE algorithm (GICE) as disclosed by Delingnon et 
a!., in "Estimation of generalized mixtures and its application in image 
segmentation" (IEEE Transactions on Image Processing, 6(10): 1364-1 375, 1997) 
can be used. GICE generally provides parameter estimates for mixtures composed 
of a various number and type of statistical distributions. 

[0214] A fifth non-restrictive illustrative embodiment of the method and 

device according to the present Invention will now be described. For the sake of 
brevity, only the differences between the fifth and first Illustrative embodiments will 
be described herelnbelow. 

[0215] In this fifth non-restrictive illustrative embodiment, the 

segmentation fast marching method allows to treat and analyze, In addition to the 
volumic information obtained from the boundary layer detections of the blood 
vessels, dynamic data retrieved from IVUS pullbacks defining a fourth dimension. 
The dynamic data generally relates to the cyclic pulsation occurring in the blood 
vessels. 

[0216] The cyclic variations of the vessel dimensions combined to 

cardiac motion (for coronary IVUS) was described, in the literature, as the 
sawtooth artifact which is generally visible on longitudinal view of the IVUS volume 
and generally caused by the blood vessel pulsations. As shown by the arrow in 
Rgure 16, the sawtooth artifact is also present in IVUS data of femoral arteries 
even without cardiac motion. In common femoral arteries for Instance, diameter 
measurements generally vary from 6.8 mm in diastole to 7.2 mm in systole for 
patients with lower limb peripheral vascular disease, as disclosed by Tai et al. in 
"In vivo femoropopliteal arterial wall compliance In subjects with and without lower 
limb vascular disease." (Joumal of Vascular Surgery, 30(5):936-945, 1999). 
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[0217] Electrocardiogram-gating (ECG-gating) acquisition was 

proposed by von Birgelen et al,, in "ECG-gated three-dimensional intravascular 
ultrasound, feasibility and reproducibility of the automated analysis of coronary 
lumen and atherosclerotic plaque dimensions in humans" (Circulation, 96{9):2944- 
2952, 1997) to remove this artifact. This is generally accomplished by acquiring 2D 
IVUS frames at a precise moment of the cardiac cycle, commonly at the end of 
diastole, which generally gives more accurate and reproducible volumic 
measurements, as disclosed by von Birgelen et al. and by Bruining et aL, In "ECG- 
gated versus nongated three-dimensional intracoronary ultrasound analysis: 
implications for volumetric measurements" (Catheterization, and Cardiovascular 
Diagnosis, 43:254-260, 1998). 

[0218] Because ECG-gatIng hardware is generally not widespread, 

retrospective gating was proposed to remove the cyclic changes on non-gated 
IVUS pullback. Change tracking In semi-automatically detected lumen contour was 
first proposed by Nadkanml et al., in "Image-based retrospective cardiac gating for 
three-dimensional intravascular ultrasound Imaging" (SPIE Proceedings: Medical 
Imaging: Ultrasonic Imaging and Signal Processing, volume 4687, pages 276- 
284, 2002). 

[0219] Another method searched for cyclic variations in contour 

features calculated in a pre-processing step as disclosed by de Winter et al., in 
"Retrospective image-based gating of intracoronary ultrasound images for 
improved quantitative analysis: The intelligate method" (Catheterization and 
Cardiovascular Diagnosis, 61^-94, 2^84). The most recent retrospective gating 
proposed method Is based on variations of the images mean gray level values by 
Zhu et al., in "Retrieval of cardiac phase from IVUS sequences" (SPIE 
Proceedings: Medical Imaging: Ultrasonic Imaging and Signal Processing, volume 
5035, pages 135-146, 2003) which states that the bigger systolic lumen, that is 
hypbechoic, generally decreases the mean gray level value of the image. 
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[0220] Some measurements can be made from the cyclic vessel 

variations. It was demonstrated by Shaw et al., in "Determinants of coronary artery 
compliance in subjects with and without angiographic coronary artery disease" 
(Journal of the American College of Cardiology, 39(1 0):1 637-1 643, 2002.) that 
plaque compression is related to the vessel cross-sectional compliance. Also the 
lumen cross-sectional area (CSA) difference between systolic and diastolic 
measurements was significantly greater in yellow plaque which generally consists 
of thin, fibrous cap with lipid-rich core and inadequate collagen content), than white 
plaque which consists of thick fibrous cap or completely fibrous, as disclosed by 
Takano et a!., in "Mqchanical and structural characteristics of vulnerable plaques: 
Analysis by coronary angloscopy and intravascular ultrasound" (Journal of the 
American College of Cardiology, 38(1 ):99-1 04, 2001 ). 

[0221] Thus, the cyclic pulsation contains information about volumic 

changes of the blood vessel wall that Is generally lost when the acquisition is ECG- 
gated. The vessel pulsation information from non-gated acquisition may be kept 
and used to reconstruct the vessel wall in 3D, at different moments of the cardiac 
cycle. With this fourth-dimensional reconstruction of the vessel wall, volumic 
accuracy and reproducibility can be achieved for measure^ments made on 3D 
image sets at specific moments of the cardiac cycle. 

[0222] To perform 4D reconstruction of the blood vessel wall, detected 

boundaries from each 2D IVUS frames first have to be classified in different wall 
pulsation phases. This step may be achieved by searching periodic components in 
measurements calculated from the detected boundaries. 

[0223] Figure 17 shows the lumen area calculated from the 

segmentation results obtained from the semi-automatic PDF-based fast marching 
method as described hereinabove. A cyclic variation is visible and could be 
assessed on these retrospective measurements. However, it is also possible to 
include the wall pulsation assessment In the boundary detection process, such that 
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this information is used to help the segmentation of the IVUS pullback. 

[0224] The wall pulsation assessment may be initiated during the initial 

contour calculation procedure (Operation 231 of Figure 24), and refined when the 
boundary detection is finished. The wall pulsation is then generally divided In a 
discrete number of phases (Operation 232 of Figure 24), such as for example 
systole, beginning of diastole and end of diastole, and a label Is assigned to each 
wall pulsation phase (Operation 233 of Figure 24). The pullback 2D IVUS frames 
are then classified and assigned to the corresponding pulsation phase label 
(Operation 234 of Figure 24). 

[0225] As shown in Figure 17, the lumen area contains the cyclic 

pulsation information. Further, the area variation between adjacent frames is 
generally used to define the classification. Since the pulsation is usually periodic, 
the cyclic pulsation infonmation is used to deform the initial subset template 158 
according to variations of the expected pulsation in the initial contour calculation of 
individual IVUS images (Operation 215 of Figure 22). The assignment of the 
pulsation phase labels to following frames can also take advantage of the periodic 
information. 

[0226] At the end of the initialization process, each 2D IVUS frame Is 

identified with a wall pulsation phase label. However, these labels may change 
because, at the end of the segmentation process, more accurate lumen areas are 
calculated. The initial labels are therefore adjustable according to their initial value, 
to the variations in area difference measurements and to the expected value 
according to the periodic variation (Operation 235 of Figure 24). 

[0227] With this pulsation assessment, the 4D data set are divided in 

3D data sets composed of all IVUS 2D frames associated to a specific cardiac 
phase label and corresponding to the different phases of the blood vessel's 
pulsation (Operation 236 of Figure 24). Volumic (boundary) measurements can 
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then be made on these separate 3D data sets (Operation 237 of Figure 24). 

[0228] Although the present invention has been described hereinabove 

by way of non-restrictive illustrative embodiments thereof, it can be modified at will, 
within the scope of the appended claims without departing from the spirit and 
nature of the subject invention. 



